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DPG method with optimal test functions 
for a fractional advection diffusion equation * 
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Abstract 

We develop an ultra-weak variational formulation of a fractional advection diffusion prob¬ 
lem in one space dimension and prove its well-posedness. Based on this formulation, we define 
a DPG approximation with optimal test functions and show its quasi-optimal convergence. 
Numerical experiments confirm expected convergence properties, for uniform and adaptively 
refined meshes. 
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1 Introduction 

In this paper we develop a discontinuous Petrov-Galerkin (DPG) method with optimal test 
functions for a one-dimensional fractional advection diffusion problem of the form 

—DD^~^Du + bDu + cu = f on/:=(0,1), 
u(0) = u(l) = 0. 

Here, D denotes a single spatial derivative, and for a € (1,2), represents a fractional 

integral operator of order a —2. Throughout, we assume that c € L°°([0,1]), b G (^^([0,1]), and 
c-Db/2 > 0. 

Fractional advection diffusion equations have been receiving increased attention over the past 
decade as modeling equations for physical phenomena in such areas as contaminant transport in 
ground water flow [3] , viscoelasticity |28| , turbulent flow |28[ [32] , and chaotic dynamics |3I| • As 
most models involving fractional order differential equations do not have closed form solutions 
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particular attention has been paid to the development of numerical approximation schemes for 
these equations. Two phenomena of fractional order differential equations which impact their 
numerical discretization and approximation are: (i) the fractional differential operator is nonlocal 
(leading to a dense coefficient matrix), and (ii) the (typical) low regularity of the solution (leading 
to slow convergence of the numerical solution to the true solution). 

The first approximation methods investigated for fractional order differential equations were 
hnite difference schemes proposed by Liu, Ahn and Turner [26], and Meerschaert and Tadjeran 
|29| . (see also |33l[ini|35|). Subsequently, finite element naiMiiiTiiioiisi] and spectral methods 
(ISllMlIig have been developed for the approximation of fractional order differential equations. 
We note that a finite difference approximation using the Griinwald formula on a uniform mesh 
leads to a Toeplitz like matrix which significantly reduces the storage required for the coefficient 
matrix, and whose linear system can be very efficiently solved using a fast Fourier transform |35j . 
Fractional diffusion problems are inherently difficult to analyze and with our method we open a 
way to deal with singularly perturbed cases (not considered here). In fact, principal objective of 
the DPG method is to provide robust discretizations of singularly perturbed problems like convec¬ 
tion diffusion naiziisiiii and wave problems |43) . The DPG method with optimal test functions 
has been developed by Demkowicz, Gopalakrishnan and co-workers. In its most common form it 
combines several ideas. These are ultra-weak variational formulations (cf. |14([8]) with additional 
trace and flux unknowns (cf. Hj), and the utilization of specific test functions which are designed 
for stability (cf. the SUPG method in |23j and test functions in |2|). Demkowicz and Gopalakr¬ 
ishnan combine these ideas in a discontinuous setting and by employing problem-tailored norms. 
Appropriately combined, the resulting DPG method with optimal test functions delivers robust 
error control and also gives access to localized a posteriori error estimation (or rather calcu¬ 
lation). For details we refer to [111 I12| . In this paper we follow precisely these steps to deal 
with equations involving fractional diffusion. By writing ([T]) as a first-order system, cf. (USD , we 
develop an ultra-weak variational formulation in Section 12.41 below. While a weak formulation 
of (HD leads to a non-symmetric, coercive bilinear form, for the DPG method with optimal test 
functions the resulting variational formulation is always symmetric, positive definite, implying 
existence of a unique solution. This is the central result of the DPG method with optimal test 
functions, stated below in Theorem[TJ Necessary conditions for its application are the well-known 
Babuska-Brezzi conditions ([2D, which we check in Section [3] for our ultra-weak formulation. A 
central step will be to extend Riemann-Liouville fractional integral operators to negative order 
Sobolev spaces and prove their ellipticity. To that end, we extend recent results from |24| . In 
our main result. Theorem [71 we show well-posedness of the underlying ultra-weak variational 
formulation and quasi-optimal convergence of the discrete scheme. In particular, we will gain 
access to error control and adaptivity. In Section [H we report on several numerical experiments 
that illustrate convergence orders of variants with uniform meshes and with adaptively refined 
meshes. 

We note that in m the authors propose a simplified Petrov-Galerkin method with optimal test 
functions for fractional diffusion. They stick to discrete spaces with continuous functions and 
calculate test functions globally. In contrast, we develop the fully discontinuous variant that 
allows for local calculations of test functions. This is particularly important for fractional-order 
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problems where inner products are defined by double integrals so that global calculations are 
prohibitively costly. Let us also mention that there is DPG-technology available for hypersingular 
integral equations [221 El]. Hypersingular operators are of order one with energy spaces of order 
1/2. For closed curves/surfaces, DPG theory can be established with integer-order Sobolev spaces 
and is then simpler in a certain way. For open curves/surfaces however, one has to return to 
non integer-order spaces. The case of hypersingular operators can be seen as a limit of fractional 
diffusion operators with orders between one and two, as considered in this paper. 


2 Mathematical setting and main results 

We use the widespread notation A< B to denote the fact that A < C ■ B where C > 0 does not 
depend on any quantities of interest. By A 2 :^ B we mean that both A < B and B < A hold. 
Throughout, suprema are taken over the indicated sets except 0. 


2.1 DPG method with optimal test functions 

We briefly recall the premises and results of the DPG method with optimal test functions, 
cf. [Ill [T2l 03]. Given a Banach space U, a Hilbert space V, and a bilinear form 6 : t/ x P —)■ R, 
we consider the following three conditions: 


b{u, v) = 0 for all r) € P 
there is a positive constant Gjnfsup such that 


rr = 0; 


^ II II ^ 6(u, r>) „ ,, 

Ginfsup||i’||y < sup — for all 17 G P; 


leU \m\u 


there is a positive constant Cb such that 


(2a) 


(2b) 


b{u, v) <CM {/||r7||y for all u £ U,v £ V. 


(2c) 


Define the so-called trial-to-test operator 0 : [7 —>■ P by 

{Qu,v)v = b{u,v) for all v £ V. 


(3) 


The following result is central to the DPG method and is, in the end, consequence of the Babuska- 
Brezzi theory [HEllM], cf. im and related references given in the introduction. 


Theorem 1. Suppose that (I2ap - (l2cll hold for a Banach space U, a Hilbert space V, and a bilinear 
form b : U x P —)• R. Then, an equivalent norm on U is given by 


\U\\e ■= sup 


b{u, v) 


^eT/ ll'*^l|y 


and Cinrsupllwllc/ < \\u\\e- 
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Furthermore, for any i € V, the problem 


find u such that b{u,v) = i(v) for allv^V 
has a unique solution, and 


\u\\e < 


\V' 


In addition, if C4ip C U is a finite-dimensional subspace, then the problem 

find Uhp G C4p such that b{ui,p,v^p) = i{v^p) for all ^»hp G 0(C/hp) 
has a unique solution, and 

\\u-ui,p\\e= inf \\u-uL\\e. 


(4) 

(5) 

( 6 ) 

(7) 


2.2 Sobolev spaces 

For s G M with s > 0 and an open interval M = (a, 6) C R, the Sobolev spaces H^{M) are 
defined via distributional derivatives and the Sobolev-Slobodeckij seminorm | • \h‘>{m) norm 
II ■ ll_ff'>(M)- Ths space H^{M) is defined as the space of functions whose extension by zero is 
in Ff*(R). The space denotes the topological dual space of while 

denotes the dual of For a finite partition T of / = (0,1) into open, disjoint, and 

connected sets, we define H^{T) := OtgT77^(7")) or, likewise, := OtcT77^(7")) with 

product norms := ll^l'rllH»(T)' write H~^{T) or H~^{T) for the duals 

of product spaces. By N := ffT we denote the number of elements in the partition and for 
V G 1/2 < s, we define the jump [u] G R-^^^ as the vector of the differences of the traces 

of V on the elements to the right and to the left of all nodes x = T_ fl T+. For the boundary 
nodes (i.e., 0 and 1), we just take traces. For v G 1/2 < s, we also define the average 

{u} as the vector of mean values of the traces of v on the elements to the right and to the left 
of all nodes. We will need certain results for this kind of spaces. From now on, we assume that 
partitions are quasi-uniform, i.e., for all T G T holds |T| ~ N~^ for N := fi^T being the number 
of elements in the partition T, and the constant involved is independent of T. We denote by 
the T-piecewise distributional derivative. 

Lemma 2. The following statements hold with constants which only depend on s: 

• Let s G (0,1/2). There holds 

\\v\\h‘>(I) ^ for all v G (8) 

• Let s G (1/2,1). There holds 

\\Drv\\^s-r^j) < N^-^\v\es^t) for all v G H^iT). (9) 
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for all V G W{T). 


( 10 ) 


• Let s G (1/2,1]. There holds 

IHI ^ 


Proof. The first statement is seen as follows: First, for T a reference interval with fixed diameter, 
there is a constant Cs > 0 such that cf. |18| and |19l Proof of Lemma 5]. 

Second, scaling arguments show that ^ T G T. Now, 

ll-llLm S llolls.,/) £ E ll>’lls.,T). (11) 

TgT 

where the second estimate follows from m Lemma 20]. To show the second statement, we 
proceed as before and use an affine transformation on every element T G T, 

ll^^ll < N\\Dv\\ 1^_,^^^ < N\\Dv\\l^_,^^y 

Here the second estimate follows as the norms involved are dual to norms on which we can 
use |18| and |19[ Proof of Lemma 5]. A quotient space argument on the reference element T, 
cf. |20j . shows 

~ 11 ^ + 

The second statement follows by application of the scaling argument 
The third statement follows easily from, e.g., 

\v{x)\ < 

Here, for example, x = T_ nr+, the second estimate follows by the Sobolev Embedding theorem, 
and the third one again by a scaling argument. □ 


Lemma 3. There holds 


Ml 2 {I) ^ \\DrT\\L 2 {i) + for all r G H^{T), 

and the hidden constant is independent ofT. 

Proof. Let (f G H^{I) be the weak solution of —D'^f = r. Then Dcf G H^{I) with distributional 
derivative = —r, and integration by parts yields 

(T,r) = -{T,D^(f)) = {DrT,D(j)) + ([r] ,D(/). 


Cauchy-Schwarz and Lemma [2l eq. (|10p imply 

By construction, < 11 t]|^ 2(/)) which concludes the proof. 


□ 
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We will need the following result on fractional seminorms. 
Lemma 4. Let s G (0,1) be fixed. There holds 


\u\h‘>{i) ^ \\Du\\hs-^i) for all u G H^I). 

where Du is the distributional derivative of u. The hidden constant does not depend on I. 

Proof. As u G ^ 2 ( 1 ), it holds Du G We can write u = Dfi + c with c G M and if G 

where ||' 0 ||/^i( 7 ) ^ definition of the distributional derivative we see 

\{u,Dif)\ = \{Du,iIj)\ < \\Du\\h-i(i}U\\h^j) < \\Du\\H-i(i)\\u\\L2ii) 

We conclude that for u G ^ 2 ( 1 ), it holds 

ll^llLf/) = (u^Dif) + (u,c) < \\Du\\h-i(i)\\u\\l^(^i) + (u,c). 

Now we apply this estimate to u — tt, where u denotes the mean value of u, and obtain 

h-u\\L2{I) ^\\Du\\h-1(I). (12) 

The standard Poincare inequality states that 


\\u-u\\Hi{i)<\\Du\\L^(^iy (13) 

The H^{I) norm can equivalently be obtained by the K-method of interpolation, cf. |34) . via 


\u - ^ ||u - u\\‘^ 


nL2{l),HfiI)]s,2 
Using (fT^ and (fTdIi . we obtain 


f 


i-2s 


inf , \\u-u- + t\\v\\Hi(i) 

v^H ^{1} 


dt 


inf 


U-U-v\\L2iI)+t\\v\\Hl(I) 


< 

inf 


veHfii) 


v=0 

< 

inf 


veHfii) 


v=0 


\\U - U - v\\l 2{I) + t\\v\\Hl(I) 
\\Du - Dv\\h-i^i) + t\\Dv\\L 2 ^j) 


Next we use that for w G ^ 2 ( 1 ) there is a ip £ H^{I) with ip = 0 such that Dip = w. We conclude 

H-pi) + t\\'w\\L2(i) 

By definition, the right-hand side is 2 ’ i® equivalent to \\Du\\‘^s-i(^jy 

This concludes the proof for a specific I. The hidden constant does not depend on I, which can 
be shown by scaling arguments. □ 


t^I |2 


\u — u 


H»(I) 


< 


f 


.-2s 


inf \\Du — w\ 
w&L2{I) 
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2.3 Fractional integral operators 

The fractional integral operators that we will use are of so-called Riemann-Liouville type. For 
/3 > 0 we denote by and the left and right-sided versions of these operators, defined 

on / = (0,1) by 

qD~^u{x) := j {x - s)^~^u{s) ds and D^'^uix) := j ~ 

We also abbreviate D~^ := qD~^. A standard textbook on this kind of operators is m- 
Recently, classical results regarding boundedness and ellipticity of these operators were extended 
in |15[ 124) . In order to obtain a variational formulation suited for DPG analysis, we need to 
extend these operators to negative order Sobolev spaces and show their ellipticity. To this end, 
let T” denote the Fourier transforms on the space 5'(M) of tempered distributions, cf. |30l Chapter 
7], defined by 

Tu{^) := (27r)“^/^ f u{x)e~^^^ dx. 

Jr 

Choosing a space of test functions which is invariant under the action of D~^ and D~^, these 
operators can be extended to the associated spaces of distributions, cf. |31) §8]. In the present 
setting, a different argument can be used. 

Lemma 5. For every s G R with —/? < s and /? > 0, the operator D~^ can he extended to a 
bounded linear operator D~^ : H^{I) 

Proof. For 0 < s, the statement was shown for and in Theorem 3.1 of |24| . It therefore 

remains to consider —/3 < s < 0. We will show the statement for s = —/3, the remaining cases 
follow by interpolation. We already know that : ^ 2 ( 1 ) —>■ H^{I) is a linear and bounded 

operator. According to ED Corollary of Thm. 3.5], it holds that 

{D~^u,v) = {u,Df^v) for all u,u G L 2 (/). (14) 

Hence, the right-hand side of ffT4)l extends D~^ to a linear, bounded operator D~^ : H~^ —> 

L2(/). □ 

Lemma 6. The operator D~^ is elliptic on 77“^/^ (T) for 0 < f3 < 1. 

Proof. For a test function ip G 77(1) holds cf. ED Thm. 7.1]. 

Then, a short computation (cf. [151 Proof of Lemma 2.4]) shows 

+ i sin(-7r/3/2) J {if)~^^‘^F{p){iff)-h/‘^F{(p)di 
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As the left-hand side of this identity is real, the imaginary part on the right-hand side vanishes. 
Furthermore, cos(—7r/3/2) > 0 for 0 < /3 < 1. We obtain 




L2{ 


The right-hand side is equivalent to the norm 
ellipticity on Since on it holds 

cf. (|lip . the proof is finished. 


h-i^/^(r)- ^ density argument shows the 


\H- 


> 

rN_> 




□ 


2.4 Ultra-weak formulation and main result 


We write m as first-order system 


a — Du = 0, 

9 (15) 

-DD'^-^a + bDu + cu = f. 

Then, we multiply these equations with r respectively v, integrate by parts piecewise on a 
partition T and rename the appearing boundary terms of D^~^a and n by a and u to obtain 


{a ,t) + {u,DrT) - {u,[t]) (16a) 

, Dtv) + {ba , v) + {cu ,v) - {a ,[v]) = {f ,v). (16b) 

The left and right-hand sides of the preceding equations dehne our bilinear form and linear form 
via 

b{u, v) := b{a, u, a,u]T,v) := {a D-j-v -|- bv) + {u , D'fT + cv) — (tt, [r]) — (a , [u]), 

£{v) := i{T,v) := {f ,v). 

Here and from now on, ^ Ff^~“/^(/) denotes the conjugate of D^~^. Define 

Ua := X L 2 {I) X x and U := H\T) x where N is the number 

of elements of T, with product norms 

II^^IlL — lkll|o/ 2 -i(j) + II^IILh) + lV“^(|ap + |up), and 

W'^Wvc ■= ll'^ll^i(r) + ll^llH“/2(r)- 

By I • I, we mean the usual Euclidean norm. Our ultra-weak formulation now reads as follows: 
given G V)(, we aim to find u ^Ua such that 

b{u,v) = £(v) for all v € Va- (17) 

For a discrete subspace t/hp C Ua, the DPG method with optimal test functions is to hnd 
iihp G Uhp such that 

5('Uijp, Uhp) — £{v\ip') for all u^p G ©^(tdip), (18) 

where ©„ : C/q —?> U is the trial-to-test operator associated with 6, cf. ([3]) . The following theorem 
is the main result of this work. It states unique solvability and stability of the continuous and 
discrete formulations (fT7|) and (IlSp . as well as a best approximation result. 


Theorem 7. For a G (1)2), / G ^ 2 ( 1 ), and arbitrary partition T, the variational formula¬ 
tion m has a unique solution u ^Ua, and 


I'^Wuc ^ II/IIl2(/)- 


Furthermore, the discrete problem (HHD has a unique solution Uhp G Hhp; and 


\u 


-«hp||f/„ < , , inf^, 

('^hp ’“hp ’ '^hp ’“hp ) ^ ^hi 


“/^lk-C^hp|li?«/2-l(j) + \\u-Uhp\\L2{I)) 


Proof. We are going to apply Theorem [T] hence we check (I2ap - (f2cl) . The condition (|2ap follows 
from Lemma [9l The condition (|2^ follows from Lemma |8l It remains to check condition 
To that end, observe first that 


sup = (ik + + ?)u||^i-o/ 2 (^) + WDtt + cvWl^^j) + N^{\[t]\ + |H|)2) 


u&Uc \m\u, 


1/2 

(19) 


For given v = {t,v) G Va we define ti G H^{I) and vi G Lf"/^(/) as the solution of Lemma [lo] 
with data F := D-pr + cv and G = r + Dj-v + bv, and write t = tq + ri and v = vq + vi. 

The functions tq and vq then fulfill the assumptions of Lemma [TTJ The triangle inequality and 
Lemmas m and m show 

ll^lltL ^ Ik + D^'^-^>Drv + + \\DrT + cv\\l^(i) + iV^/^(|[r]| + |[u]|). (20) 

The equations m and poll show condition (I2bp . Theorem [T] shows that there are unique 
solutions u and Uhp of the problems (I17h and fllSp which fulfill stability ([5)) and best approxi¬ 
mation and that 


C'infsupllwllf/^ < llwlbc 


sup 

VGVa 


b{u, v) 

Ikllv 

M lira; 


Lemma [8] shows that 


,inf \\u - u'^pWe^ < 

r ^ - < p \\ hc >/ 2 -^ j ) + \\U - <pI|l2(/)) ■ 

«p)<p)<p)<p)et^hp ^ ^ 

Here, owing to the fact that u and a are just finite-dimensional vectors, the functions and 
can be omitted on the right-hand side. □ 
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3 Technical results 

The first lemma states boundedness of the bilinear form b. 

Lemma 8. For a € (1, 2), 

1 /2 

\b{u,v)\< + l|^;||y„, 

with a constant independent ofT- In particular, \b{u,v)\ < Cb||u||;7„||r’||v'a) where the constant 
Cb depends on N. 

Proof. By Lemma [2l eq. (|10p . we have 

and \{a ,M)\ ^ 

The triangle inequality, Lemmas [2] and [5l the definition of ^ g L°°([0,1]), b € ^^([0,1]), 

and 1 < q: show 

||r + < ||r||j^i-a/2(/) + H-Dt-uH j^c,/2-i(/) + \\bv\\fji-a/2(^i^ 

< (ikllnTT) + ll^llH“/2(r)) 

and 


WDtt + cu||i2(7) < llrll^qr) + lbll//-/2(r)- 

We finish the proof with the triangle and Cauchy-Schwarz inequalities. □ 

Lemma 9. It holds that 


b{u, v) = 0 for all V ii = 0. 

Proof. The direction is clear, and we proceed with the implication =^. Using r € C^(/) 
in ()16ap shows that the distributional derivative of u fulfills Du = a. As cr G H°‘/‘^~^{I), we 
conclude that u G In a second step, using functions r G C°°{T) for all T G T in (I16al) 

and integrating by parts shows that ft = rt at inner nodes as well as u{a) = u{b) = 0. Hence 
u G We plug in fj = Du in ()16bp and obtain the variational formulation 

{D°‘~'^Du , Dv) + {bDu , v) + {cu ,v) = 0 for all v G H°'/‘^{I). 

According to m Section 3], the bilinear form on the left-hand side of this formulation is elliptic 
on H°^/‘^{I). We conclude that u = 0 and hence <7 = 0. Then, u = 0 and a = 0 follow 
immediately. □ 
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3.1 Analysis of the adjoint problem 

Lemma 10. For F € ^ 2 ( 1 ) and G € there exists a solution r € v G 

of 


Dt + cv = F 
T + D^oc-‘^'>*Dv + bv = G 


(21) 


such that 


ll^ll//“/2(/) + lkllHl(/) ^ l|-^llL2(-f) + I|G'||_H-1-“/2(J-)- (22) 

Proof. Consider the variational formulation to find v G such that 

{Dv , D^-^Dcf) + {bv , D(j)) + {cv ,(j)) = {F,(j))- {DG , (f) for all 0 G 

According to m Section. 3], the bilinear form of this formulation is elliptic on The lin¬ 

ear functional on the right-hand side is bounded in if"/^(/) with constant ||T||£, 2 (/) + ||G||^i-c</ 2 (/p 

Hence, there exists a unique solution v G i/“/^(/) which satisfies 

\\ v \\ hc ,/ 2 i ^ j ) < II-^IIlzH) + I|G'IIhi-“/2(/) 

Now define r := —D^°‘~^^*Dv — bv + G. A priori, r G but the definition of v shows 

{T,D(j)) = {cv-F,cl)) for all (/. G C^(/). 

Hence, r G H^{I) and Dt = F — cv. The bounds on r follow immediately. □ 

Lemma 11. Suppose that r G H^{T) and v G fulfill 

D-j-t -|- cu = 0 (23a) 

T + D^‘^-^>Drv + bv = 0, (23b) 

on /. Then 

Ikllfl'i(T) + lbll//“/2(r) ~ (INI + INI) • 

Proof. We proceed in three steps. 

Step 1 : Let if G il“/^(/) be the unique variational solution of —DD^ + bDip + cif = —v, 
cf. [m Section 3], i.e., 

, D(P) + {bDif ,(f)) + {ci;,(p) = -{v, cf) for all 0 G 
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(24) 


so that ||'0||/^a/2(7) ^ II^IIl 2 (/)' Due to Lemma [5] and a — 2 < a/2 — 1, it holds 

The equation solved by ^ implies that the distributional derivative of is given by 

DD°^~^D'ijj = bDil; + + u G such that D^~^D'ijj G Using Lemma 01 we 

see 

< \\DD^-^D'il;\\fjo./2-^i^ = WbUip + ciIj + v\\ho./2-ii^j) 

(25) 

< ||V’IIhc./2(j) + ||?;||l2(/) ^ \Hl2{i)- 

We may also integrate by parts and use (|23p to obtain 


{v ,v) = —(ZD" , D^-v) — {bDtp , v) — {ctp, v) + (ZD" , [u]) 

= {Di/j , -ZD("-2)*ZDru - 6u - r) + (V’, [r]) + {D^-^D'ijj, [u]) 
= (V'Jt]) + (D""^ZDV’, H). 

Lemma [21 eq. m , estimates (1241) . (1251) . and stability of ijj yield 

ll*'llLa)S"'''(IWI + IHI)ll»ll«"/-(T)- 

Step 2 : Piecewise integration by parts shows 

{D']-{bv), v) = (vDb , v) + [bD-yv , v) 

= {vDb,v) - {v,Dr{bv)) - (H ,{6u}) - ({u}, [6?;]), 


which gives 


(26) 


{Dr{bv) , v) = (vBb/2 , v) - l/2([u], {bv}} - l/2({u} , [bv]). (27) 

Now we multiply (|23ap with v and insert (I23bp as well as ([271) . Then, as G ZZ"/^(T) 

by (|23bp . integration by parts gives 

0 = {Drv , D^-^Drv) + {v{c - Db/2), u) 

+ ([ZD("-2)*ZDru], {u}) + ({ZD("-2)*ZDru} , [u]) + l/2([u], {bv}) + l/2({u} , [bv]) 

As c— Db/2 > 0 and ZD"“^ is elliptic in ZZ"/^“^(T) due to Lemma[6l we obtain with the triangle 
inequality 

l|DrHI^./2-i(r) < KM > W)l + KM, H)l + KH,{M)l + l(N, W)l- 
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All terms on the right-hand side of this inequality are treated with Lemma [2l eq. (fTOl) . For the 
second term, we additionally use Lemma [3] and (I23all and get 


KW.H>l S iv‘''"l|TllH.(T) ■ IHI < iv'''"l|!>llL.(;)IHI + iviMI ■ IHI 

< ■ IWI- 

We conclude that 


|i^lH«/2(r) - \\DTv\\\c/2-i^r) ~ (I Ml + I Ml) MIM“/2(r)> (28) 

where we have used Lemma 0] for the first estimate. Adding (1261) and (|28p and dividing by 
\\v\\ho./ 2 ^j^ gives 

llollH.flfnSJvs/qilTlI + lHI). (29) 

Step 3: It remains to show the bound for r. As r € ^ 2 ( 1 ), we can write r = + t with 

ijj G H^{I) and t G M such that ||'i/’||//i(/) + |^| ^ II'^IIl 2 (/)' Integration by parts, identities (I23p . 
Lemma [2] eq. cni), and Cauchy-Schwarz show 


(r, r) = (cv , V') + ([r], V’) - '^'>*Drv + bv,t) 

< (llKlIfaU) + Jv‘'"I|t 1| + (D<“-2>*Drt..l)) l|T||fa,„. 

For the last term, we use = x^““/r(2 — a + 1), cf. |311 Section 2.5], and integration 

by parts to compute 






{[v]y-n 

r(2-a + l) 


(31) 


< \\v\\hc^/2^j^ + N^/\v\\. 

Here, the last estimate follows by direct computation. Combining the estimates (1301) . (I3ip . 
and (j29p . we obtain 


l|r|lw/)<iV“''"(l|Tl| + |MI). 

An estimate for D-yr is obtained from (|23ap and (|29l) . This concludes the proof. □ 


4 Numerical Examples 

4.1 Discretization and approximated optimal test functions 

Let us briefly hx some notation: We consider the discrete subspace 

t/hp(T) := UP{T) X U\T) X X C C/„, 
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where 


U^{'T) := {f G L 2 {I) : v\t is polynomial of degree at most p VT G T} 


is the space of T-elementwise polynomials of degree p G Nq. Note that dim([/hp(T)) = (p + 
q + 4)A^. Given a basis {uj | j = 1,..., dim(f7hp(T))} of f7hp(7”), the optimal test functions 
&{uj) G (j = 1) ■ ■ ■) dim(C/hp(T)) are computed by solving the problems 

{Q{uj),v)v^ = b{uj,v) for all v € Va = H\T) x (32) 


For V = {t,v),w = {p,w) G Va the f^-inner product is given by 

, V / ^ ^ N / N f f (v(x) — v(y))(w(x) — wiv)) 

{v ,w)v„ = + {DrT,Drp)i + {v,w)i + ' 


tgT 


dy dx, 


T JT 


X - 

on Va- Since the dehnition 
of the optimal test functions (f32ll involves the inhnite-dimensional space Va, we approximate 
©^(wj) G Va by Qa,h{xLj) G 14p(7”) := U'^{T) x U"'{T) with m,n € Nq, i.e., instead of (1321) we 
solve for j = 1,..., dim(f7hp(7”)) the problem 


which induces our chosen local norm = ll'^ll^i( 7 ') + ll^ll^a/ 2 (- 7 -) 


{V)a,h{uj) ,Vk)vc =b{uj,Vk) A; = 1,... ,dim(14p(T)). 


(33) 


The inner product {v ,w)y^ is computed analytically for functions v,w G Vhp(T). It is seen 
immediately that choosing m and n too small in comparison with p and q leads to a system 
which is not well posed. This question is investigated in HZ]. The authors show that in the case 
of the Poisson equation in and p = q, using polynomial degrees n = m which are higher than 
p + d is sufhcient in order to obtain well-posedness and best approximation results. 

Altogether we have to assemble the matrices B := (B^j) and 0 := {@ki) with 

Bkj := b{uj,Vk) and @ke ■= {ve ,Vk)Vc, 

where Uj and Vk, j = 1,... , dim(I7hp(T)), k = 1,... , dim(Idip(T)), are the basis functions 
described above. Note that 0 has a sparse structure, whereas B contains a dense block cor¬ 
responding to the discretization of the fractional integral operator. With the dehnition of the 
right-hand side vector 


tj := i{vj) for all j = 1,..., dim(I4p(T)) 

the computation of the DPG solution ([6|) consists in solving the linear system 

B^0“^Bx = B^0~^f. (34) 

An advantage of the DPG method is that, by design, we can evaluate the error in the energy 
norm. We dehne the local contributions of the error in the energy norm on an element T G T, 
est(T), as 

est(r)2 := ^ (f - Bx), (0“i(f - Bx)) . (35) 

{j:Vj\T/=0 for T'j^T} 


14 



Then, with rh G Vhp denoting the element corresponding to the vector 0 ^ (f — Bx) it holds 


est^ := ^ est{Tf = (f - Bx)^ (0~^(f - Bx)) = 
Ter 


2 

Vh- 


(36) 


Let us discuss the convergence rates we can expect. Due to standard approximation results of 
the L 2 -orthogonal projection TTp : L 2 {I) U^{T) we have 




^ 11^ ^ ^11 ^ AT—min(p+l,s) I 

< ||(T - 7rp(T||i2(7) < A/ ’ ^1 




and 


inf 

“hpGf^nr) 




According to Theorem [71 this yields 

11^ - «hp||f/„ < est < iV--n(,+i,p+i,.,.) . ( 37 ) 

Here, the fact that est can be included in this estimate in this way follows from Theorem [H For 
the numerical examples where the exact solution u = {a, u, a, u) is known, we can compute the 
exact error ||w||Lfa- For this we define the quantities 


err(u/i) := ||u - Uh\\L2{l), 

err{ah) := - ah)\\L2{i), 

err({t/i) := N~^^'^\u - Uh\, 

err(5fe) := - ah\. 

Here, u are the evaluations of the function u at the interior nodes (i.e. without the endpoints of 
the interval I = (0,1)) of the mesh T and a are the evaluations of D^~°‘Du at the nodes of T. 
We emphasize that the norms err(fT/i), evr{uh), and (ah) to measure the error of approximations 
to u, u, and a are stronger than those contained in the norm ||rr||(7Q on the left-hand side of (I37p . 
However, the experiments show that we have optimal convergence rates also in these stronger 
norms. We emphasize that we even have the rigorous error bound 

11“ “ “hpIlL - est^ - err(fj/j)2 erv{uhf. 


4.2 Example 1 


We consider the following example, see also m Section 5, Example 2]: Let I = (0,1), a = 2)12, 
b{x) := 1/2, c(x) := 1/2 for x G I and prescribe the exact solution u{x) = — x^. Then, the 

right-hand side is given by 


fix) = -2 


r( 2 ) 


„2-a 


r(3) 


r(3 — a) 


r(4-a) 


— -x^ — X^ -|- X. 
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Furthermore, straightforward calculations show 

a{x) = Du = 2x — 3x^, 

D^-^a{x) = D^-^Du = 2—- 3— 

r(4 — a) r(5 — a) 

We consider uniform meshes on I with mesh-size h = 1/N and N = ^T. Figure [1] shows 
results for different values of p, g, m, n. As u and cr are both smooth, we expect from (|371l 
that est = and the numerical experiments reflect this expectation. We 

even see in the experiments the simultaneous approximation orders err(tt/i) = and 

err(cj/i) = The trace errors err(u/i) and err(CT/i) show higher convergence rates in 

all cases. In the case p = 0,q = l,m = 2,n = 2 (upper right plot), est converges slightly faster 
than err(iT/i) but slower than err(u/i). 


4.3 Example 2 

For the next example we prescribe the exact solution u{x) = x^ — x with 1/2 < A < 3/2 on 
I = (0,1), see also |15[ Section 5, Example 3]. The right-hand side as well as a are given by 


f{x) = 


+ j.A-« , ^ i-g 

F(A + l-a) r( 2 -a) 


cr(x) = Du = \x^ ^ — 1, 


D^-^Du{x) 


r(A-|-2 —a) F(3 —a) 


We have u G and a G for all e > 0, and hence, due to 1/2 < A < 3/2, 

with a view to (I37p . we expect a convergence rate of est = However, with a view to 

the norm || • ||[/^, the expected rate, dictated by a in this case, would be This 

is what we will see for uniform refinement. In order to regain the optimal convergence orders 
we utilize an adaptive strategy where we use est(T) as local rehnement 
indicators and mark elements Ad C 7~ according to Dorfler’s marking criterion 


0est^ < ^ est(T)^, 
TeM 


(38) 


where we use 9 = 0.4 and A4 is a set of minimal cardinality. Note that 9 = 1 means uniform 
refinement, i.e., A4 = T. Each marked element T G A4 is bisected such that local quasi¬ 
uniformity 

diam(T) 

max -—— < 2 

T,T'gr diam(T') 

TnT^7^0 

is preserved. Figure [2] shows est and the error quantities for the parameters 


A = 0.6, a = 1.2. 
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In the upper left plot the results for uniform refinement and p = q = 0, n = m = 2 are given. 
We observe the convergence rate est = expected, also for 

the separated error contributions, we observe reduced convergence rates. Adaptive refinement 
recovers the optimal rate as is seen in the three remaining plots. As in 

Example 14. 2 [ we see that the traces even have better convergence rates. 

4.4 Example 3 

In the last experiment we set f{x) := log(x) for x G / = (0,1) and note that / G L 2 {I). For this 
right-hand side we do not know the explicit form of the solution u. Therefore, we only plot the 
error in the energy norm est for different values of p, q, m, n and a, respectively. Throughout, we 
set p = q as well as m = n = p + 2. Figure [3] shows the error in the energy norm est for a = 1.6 
(left) and a = 1.8 (right). We compare uniform refinement (9 = 1) and adaptive refinement 
with 9 = 0.4 for p = q = 0. Moreover, we plot the results in the adaptive case with p = q = 1 
resp. p = q = 2. We observe that for adaptive refinement we obtain convergence rates p+1, i.e., 
est = whereas for uniform refinement we get only the suboptimal rate a/2 — 1/2. 
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Figure 1: Experimental convergence rates for Example 1 from Section [4.21 Uniform mesh refine¬ 
ment is used througout. 
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p = 0,q = 0,m = 2,n = 2,9 = l p = 0, q = 0, m = 2, n = 2, 6 = 0.4 





Figure 2: Experimental convergence rates for Example 2 from Section [4.31 Uniform mesh refine¬ 
ment (upper left) and adaptive mesh refinement (upper right and below). 
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error in energy norm 


CL = 1.6 OL = 1.8 



Figure 3: Experimental convergence rates for Example 3 from Section 14.41 The choice 6=1 
refers to uniform mesh refinement, while 9 = 0.4 refers to adaptive mesh refinement. 
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